#include <cmath>
#include <complex>

#include <Vector4.hpp>

#include <Spin.hpp>

class LS {
public:
  int L;
  Spin S;
};

inline double _dfact(double __x){
  if((__x < 0.00001) && (__x >= 0.0)) return 1.;
  if(__x < 0) return 0.;
  return __x*_dfact(__x - 1.);
}

inline int _factorial(int __i) {
  int f = 1;
  if((__i == 0)||(__i == 1)) f = 1;  
  else{
    while(__i > 0){
      f = f*__i;
      __i--;
    }
  }
  return f;
} 

#define MAX(x,y) (x>y ? x : y)
#define MIN(x,y) (x<y ? x : y)

// std::complex<double> _BreitWigner(const double px, const double py, const double pz, const double e,
//  				  const double mass, const double width)
//  {
//    std::complex<double> i(0.,1.);
//    double p2 = e*e - px*px - py*py - pz*pz;
//    return mass*width/(p2 - mass*mass + i*mass*width);
// }


double _BlattWeisskopf(Spin L,double p)
{
  const double Pr = 0.1973; // Gev/c corresponds to 1 fermi
    
  double ret;
  double z = (p/Pr) * (p/Pr);

  if (L.Denominator() != 1) {
    std::stringstream ss; 
    ss << "Nature does not allow orbital angular momentum to come in half-integer increments (`" << L.ToString() << "')";
    throw std::domain_error(ss.str().c_str());
  }

  int l = L.Numerator();

  switch (l) {
  case 0:
    ret = 1.0;
    break;
  case 1:
    ret = sqrt((2.0 * z)/(z + 1));
    break;
  case 2:
    ret = sqrt((13.0 * z * z)/(pow(z - 3.0,2.0) + 9.0 * z));
    break;
  case 3:
    ret = sqrt( (277.0 * pow(z,3.0))/(z * pow(z - 15.0,2.0) + 9.0 * pow(2.0 * z - 5.0,2.0)));
    break;
  case 4:
    ret = sqrt( (12746.0 * pow(z,4.0))/( pow(z * z - 45.0 * z + 105.0,2.0) + 25.0 * z * pow(2.0 * z - 21.0,2.0)));
    break;
  case 5:
    ret =  sqrt(z*z*z*z*z/(893025.0 +99225.0*z +6300.0*z*z +315.0*z*z*z +15.0*z*z*z*z +z*z*z*z*z));
    break;
  default:
    std::stringstream ss;
    ss << "Blatt-Weisskopf not defined for L == " << l;
    throw std::domain_error(ss.str().c_str());
  }
  return(ret);
}

double _lambda(double a, double b, double c)
{
  return( a*a + b*b + c*c - 2.0*(a*b + b*c + c*a) );
}

double _BreakupMomentum(Vector4& p, Vector4& child1_p, Vector4& child2_p) {
  double Msq  = p.M2();
  double m1sq = child1_p.M2();
  double m2sq = child2_p.M2();
  
  double lam = _lambda(Msq,m1sq,m2sq);
  return sqrt(fabs(lam/(4*Msq)));
}

double _BreakupMomentum(const double pdg_mass, Vector4& child1_p, Vector4& child2_p) {
  double Msq  = pow(pdg_mass,2);
  double m1sq = child1_p.M2();
  double m2sq = child2_p.M2();
  
  double lam = _lambda(Msq,m1sq,m2sq);
  return sqrt(fabs(lam/(4*Msq)));
}

double _BreakupMomentum(const double pdg_mass, const double m1, const double m2) {
  double Msq  = pow(pdg_mass,2);
  double m1sq = pow(m1, 2);
  double m2sq = pow(m2, 2);
  
  double lam = _lambda(Msq,m1sq,m2sq);
  return sqrt(fabs(lam/(4*Msq)));
}


std::complex<double> _BreitWigner(double m, double m1, double m2, double pdg_mass,
				  double pdg_width, int l_decay)
{
  std::complex<double> ret;

  double m0 = pdg_mass;
  double Gamma0 = pdg_width;
  int l = l_decay;
  double q      = _BreakupMomentum(m, m1, m2);
  double q0     = _BreakupMomentum(m0, m1, m2);

    
  std::complex<double> i(0, 1);

  double GammaV = Gamma0*(m0/m)*(q/q0)*(pow(_BlattWeisskopf(l,q),2)/pow(_BlattWeisskopf(l,q0),2));
  ret = (m0*Gamma0)/(m0*m0-m*m-std::complex<double>(0,1)*m0*GammaV);
  
  return (ret);
}

std::complex<double> _BreitWigner(Vector4& p, Vector4& child1_p, Vector4& child2_p,
				  const double pdg_mass, const double pdg_width, int l_decay)
{
  std::complex < double >ret;
  
  double m0     = pdg_mass;
  double Gamma0 = pdg_width;
  double q      = _BreakupMomentum(p, child1_p, child2_p);
  double q0     = _BreakupMomentum(pdg_mass, child1_p, child2_p);
  double m      = p.M();
  int l         = l_decay;
  
  double GammaV = Gamma0*(m0/m)*(q/q0)*(pow(_BlattWeisskopf(l,q),2)/pow(_BlattWeisskopf(l,q0),2));
  ret = (m0*Gamma0)/(m0*m0-m*m-std::complex<double>(0,1)*m0*GammaV);
  
  return (ret);
}

std::vector<LS> _GetValidLS(const Spin &__j,int __parity,
			   const Spin &__s1,int __p1,
			   const Spin &__s2,int __p2) 
{
  std::vector<LS> valid_ls;
  LS ls;
  
  for(Spin S = fabs(__s1 - __s2); S <= (__s1 + __s2); S++){
    for(int L = (int)fabs(__j - S); L <= (int)(__j + S); L++){
      // check the parity
      if(fabs(__p1*__p2*pow(-1.,L) - __parity) < 1.e-5) {
	ls.L = L;
	ls.S = S;
	valid_ls.push_back(ls);
      }
    }
  }
  return valid_ls;
}

double _Clebsch(const Spin &__j1,const Spin &__m1,const Spin &__j2,
		const Spin &__m2,const Spin &__J,const Spin &__M)
{

  if((__m1 + __m2) != __M) return 0.;
  if(abs(__m1) > __j1) return 0;
  if(abs(__m2) > __j2) return 0;

  // convert to pure integers (each 2*spin)
  int j1 = (int)(2.*__j1);
  int m1 = (int)(2.*__m1);
  int j2 = (int)(2.*__j2);
  int m2 = (int)(2.*__m2);
  int J = (int)(2.*__J);
  int M = (int)(2.*__M);

  double n0,n1,n2,n3,n4,n5,d0,d1,d2,d3,d4,A,exp;
  int nu = 0;
  
  double sum = 0;
  while(((d3=(j1-j2-M)/2+nu) < 0)||((n2=(j1-m1)/2+nu) < 0 )) { nu++;}
  while (((d1=(J-j1+j2)/2-nu) >= 0) && ((d2=(J+M)/2-nu) >= 0) 
	 &&((n1=(j2+J+m1)/2-nu) >= 0 )){
    d3=((j1-j2-M)/2+nu);
    n2=((j1-m1)/2+nu);
    d0=_dfact((double) nu);
    exp=nu+(j2+m2)/2;
    n0 = (double) pow(-1.,exp);
    sum += ((n0*_dfact(n1)*_dfact(n2))/(d0*_dfact(d1)*_dfact(d2)*_dfact(d3)));
    nu++;
  }

  if (sum == 0) return 0;

  n0 = J+1;
  n1 = _dfact((double) (J+j1-j2)/2);
  n2 = _dfact((double) (J-j1+j2)/2);
  n3 = _dfact((double) (j1+j2-J)/2);
  n4 = _dfact((double) (J+M)/2);
  n5 = _dfact((J-M)/2);
  
  d0 = _dfact((double) (j1+j2+J)/2+1);
  d1 = _dfact((double) (j1-m1)/2);
  d2 = _dfact((double) (j1+m1)/2);
  d3 = _dfact((double) (j2-m2)/2);
  d4 = _dfact((double) (j2+m2)/2);
  
  A = ((double) (n0*n1*n2*n3*n4*n5))/((double) (d0*d1*d2*d3*d4));


  return sqrt(A)*sum;		
}

double _Gamma(double __z) {  
  float x = fabs(__z);
  // Coefficients for the series expansion
  double c[7] = { 2.5066282746310005, 76.18009172947146, -86.50532032941677
		  ,24.01409824083091, -1.231739572450155, 0.1208650973866179e-2
		  ,-0.5395239384953e-5};

   double y = x;
   double tmp = x + 5.5;
   tmp = (x+0.5)*log(tmp)-tmp;
   double ser = 1.000000000190015;
   for(int i = 1; i < 7; i++) {
      y += 1;
      ser += c[i]/y;
   }
   double lnGamma = tmp+log(c[0]*ser/x); // natural log of Gamma(z)
   if(__z > 0) return exp(lnGamma);   
   else if(__z == 0.) return 0.;   
   else {
     double pi = 3.1415926535;     
     return (-1.*pi)/(exp(lnGamma)*sin(pi*x)*x);
   }   
}

double _Wigner_d(const Spin &__j,const Spin &__m,const Spin &__n,double __beta){
  int J = (int)(2.*__j);
  int M = (int)(2.*__m);
  int N = (int)(2.*__n);
  int temp_M, k, k_low, k_hi;
  double const_term = 0.0, sum_term = 0.0, d = 1.0;
  int m_p_n, j_p_m, j_p_n, j_m_m, j_m_n;
  int kmn1, kmn2, jmnk, jmk, jnk;
  double kk;

  if (J < 0) {
    std::stringstream ss;
    ss << "Illegal value " << __j << " for J (J >= 0)";
    throw std::domain_error(ss.str().c_str());
  }

  if (fabs (M) > J) {
    std::stringstream ss;
    ss << "Illegal value " << __m << " for M (abs(M) <= J, J == " << __j << ")";
    throw std::domain_error(ss.str().c_str());
  }

  if (fabs (N) > J) {
    std::stringstream ss;
    ss << "Illegal value " << __n << " for N (abs(N) <= J, J == " << __j << ")";
    throw std::domain_error(ss.str().c_str());
  }
  
  if (__beta < 0) {
    __beta = fabs (__beta);
    temp_M = M;
    M = N;
    N = temp_M;
  }

  m_p_n = (M + N) / 2;
  j_p_m = (J + M) / 2;
  j_m_m = (J - M) / 2;
  j_p_n = (J + N) / 2;
  j_m_n = (J - N) / 2;
  
  kk = (double)_factorial(j_p_m)*(double)_factorial(j_m_m)
    *(double)_factorial(j_p_n) * (double)_factorial(j_m_n) ;
  const_term = pow((-1.0),(j_p_m)) * sqrt(kk);	
  
  k_low = MAX(0, m_p_n);
  k_hi = MIN(j_p_m, j_p_n);

  for (k = k_low; k <= k_hi; k++) {
    
    kmn1 = 2 * k - (M + N) / 2;
    jmnk = J + (M + N) / 2 - 2 * k;
    jmk = (J + M) / 2 - k;
    jnk = (J + N) / 2 - k;
    kmn2 = k - (M + N) / 2;
	
    sum_term += pow ((-1.0), (k)) *
      ((pow (cos (__beta / 2.0), kmn1)) * (pow (sin (__beta / 2.0), jmnk))) /
      (_factorial (k) * _factorial (jmk) * _factorial (jnk) * _factorial (kmn2));
  }

  d = const_term * sum_term;
  return d;
}

std::complex<double> _ReggePropagator(double __t,double __s,double __a,double __b,
				     const Spin &__spin,int __sig,int __exp_fact)
{
 
  std::complex<double> prop,numerator,denominator;
  double alpha = __a*__t + __b; // trajectory
  double pi = 3.1415926535; 
  std::complex<double> i(0.,1.); 
  numerator = pow(__s,alpha - __spin)*pi*__a;
  numerator *= ((double)__sig + ((double)__exp_fact)*exp(-i*pi*alpha));
  double gamma_arg = alpha + 1. - __spin;
  if(gamma_arg < 0.) denominator = 2*pi/(fabs(gamma_arg)*_Gamma(fabs(gamma_arg)));
  else if(gamma_arg == 0.) denominator = 2*pi;
  else denominator = 2*sin(pi*gamma_arg)*_Gamma(gamma_arg);

  prop = numerator/denominator;

  return prop;
}

Spin _GetSpin(const std::string &__spin){  
  Spin j;
  std::string::size_type div_pos = __spin.find("/");
  if(div_pos != std::string::npos){
    std::string tmp;
    tmp.assign(__spin,0,div_pos);
    double numer = atof(tmp.c_str());
    tmp.assign(__spin,div_pos + 1,__spin.size() - div_pos - 1);
    double denom = atof(tmp.c_str());
    j = numer/denom;
  }
  else j = atof(__spin.c_str());
  
  return j;
}

inline std::complex<double> _Wigner_D(double __alpha, double __beta, double __gamma,
				      const Spin &__j, const Spin &__m, const Spin &__n){
  std::complex<double> i(0.,1.);
  return exp(-i*(__m*__alpha + __n*__gamma))*_Wigner_d(__j,__m,__n,__beta);
}

